Gradient corrections to the local density approximation 
for trapped superfluid Fermi gases 
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Two species superfluid Fermi gas is investigated on the BCS side up to the Feshbach resonance. 
Using the Greens's function technique gradient corrections are calculated to the generalized Thomas- 
Fermi theory including Cooper pairing. Their relative magnitude is found to be measured by the 
small parameter (d/ Rtf) 4 , where d is the oscillator length of the trap potential and Rtf is the 
radial extension of the density n in the Thomas-Fermi approximation. In particular at the Feshbach 
resonance the universal corrections to the local density approximation are calculated and a universal 
prefactor Kw = 7/27 is derived for the von Weizsacker type correction Kw(fi 2 /2m)(V 2 n 1//2 /n 1 ' 2 ). 



PACS numbers: 31.15.xg,74.20.Fg,67.85.Lm 



I. INTRODUCTION 

Fermi gases below the degeneracy temperature have 
been the subject of intensive research in the last years 
both experimentally and theoretically (See for reviews 
P, [![)■ Particular interest has been devoted to the pos- 
sible superfluid state whose creation and properties have 
been studied for both negative and positive values of s- 
wave scattering lengths a, characterizing the interaction 
between the particles. At the Feshbach resonance 0-0| 
a becomes infinity and certain universal behavior shows 
up. An important aspect of the problem is that the gas is 
trapped and thereby is inhomogeneous. When the energy 
gap function exceeds the level spacing near the Fermi sea 
a local density approximation (LDA) is applicable. As a 
simplest approach in its spirit neglecting the space gra- 
dients of the density and the gap function the Thomas- 
Fermi theory was generalized to include superfluid pair 
correlation results |8j , when the system is treated in the 
generalized Hartree-Fock method |9| . Since the Thomas- 
Fermi approximation is widely used in case of trapped 
gases it is desirable to investigate systematically the cor- 
rections to it, even if they are expected to be small for 
large particle numbers, except in the surface region (Here 
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the gradient corrections make explicitely visible the lim- 
its of the usual LDA results). For particle numbers, how- 
ever, which are treated in Monte-Carlo simulations the 
Thomas-Fermi theory needs corrections. More impor- 
tantly it makes possible to extend the concept of univer- 
sality at unitarity 10]. In particlular we derive in this 
paper a universal prefactor for the von Weizsacker type 
correction to the generalized Thomas-Fermi theory (see 
for a review of the von Weizsacker correction in normal 
systems [U). 

In the present paper gradient corrections are calcu- 
lated up to second order at zero temperature. Baranov 
12] studied the gradient corrections even at finite tem- 
peratures in cases when Eilenberger's equations fl3j are 
applicable. That approach is different from ours, which 
is free from this restriction. The applied technique here is 
based upon the equation of motion as expressed in terms 
of the Green's functions. The method has been developed 
first to the electron gas of the atoms, which is of course a 
normal system (T3 | . It has been generalized to superfluid 
state somewhat later independently for superconductors 
in slowly varying magnetic field 15 1 and for nuclei 
The latter work is most closely related to the present 
one. The resulting expressions are rather cumbersome, 
but considerably simplify at unitarity. To evaluate them 
we choose the mean-field BCS (MF-BCS) model intro- 
duced by Leggett, Eagles, Nozieres and Schmitt-Rink 
[16l4l8l |. which neglects the self-consistent Hartree-type 
terms. We start however, from the generalized Hartree- 
Fock (GHF) model Q to present the results in a more 
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complete form for future use. The Hamiltonian is 
+ / rf 3 ^ 3 r'^(r)^(r')«(r,r')Va'(r')^(r), (1) 

a. a' 

where U ex t(r) is the trapping potential, /i is the chem- 
ical potential, v(r — r') describes the interaction and a 
stands for the internal degrees of freedoms. We assume 
two equally populated hyperfme states and a =t, \ will 
be termed as spin. In GHF approximation the Hamilto- 
nian simplifies to 

ffmf = Z)/ ^» (-£ y2 + ^- t (r) - A*) ^a(r) 
+ ]T /d 3 r dV «(r, r')/i CT><T (r, r)V>+, (r')Vv (r') 

(T,<T' ^ 

- f d3r dV v ( r > ( r > OV^WW (r') 

+ \Hj^ T dV w ( r ' r ') (W' ( r > r ')^(r)^ (r') + He) . 



where 



(2) 



Here 



n(r)tot = Z nCT ( r ) = 2n ( r ); 

17 

ra ff (r) = /z CTjCT (r,r), 
/v, ff (r',r) = (</,+ (r)^(r')). 



(3) 



The first line in ([2]) contains the one-particle term of (JTJ , 
the second line is the Hartree term, the third is the Fock 
term and furthermore the Cooper pairing is represented 
by the last line, where 



X ^ CT (r',r) = (Vv(r)Vv(r')) 



(4) 



The correlation functions x an d h have to be determined 
self-consistently. 

We shall consider the special case when the interaction 
can be approximated by a contact potential 



v(r, r') 



8{T-v') = g5{v-v'). (5) 



In case of contact interaction the first three lines of the 
Hamiltonian ([2} can safely joined together as follows 

H mf = Yj d 3 r^(r) ("|^ 2 + U(r) - fij </v(r) 
+\llj d * r dV w ( r ' r ') (W ( r ' r ')^(r)^ (r') + He) , 

a. a' 

(6) 



C/(r) = U ext (r) +gn(r), 



(7) 



but we keep the fourth line of ((2]) as it is, because 
Xo.o' (r, r) is not a well-defined object. 

The paper is organized as follows. In Section [II] we 
present the equations for the Green's functions, while 
in Section IIIII their perturbation series are presented. 
The self-consistent scheme for the density and the gap 
function is worked out in Section IPV1 up to second order 
in h to the local density approximation, which can be 
regarded as a generalized Thomas-Fermi theory. In the 
second part of the paper the MF-BCS model is applied. 
In Section [V] the second order corrections are evaluated 
perturbatively in the case of a general external potential. 
Section IVIl is devoted to the problems of the unitary 
gas, in particular the prefactor of the von Weizsacker 
type correction is calculated. In Section IVII1 the trap 
potential is assumed to be an isotropic harmonic one to 
make some features more visible. Section fVIIII contains 
the discussion of the results. 



II. FORMULATION 

The gradient expansion can be best derived using the 
one particle normal 

G>,<r'(ri,ti;r 2 ,t 2 ) = -i(r^(n,*i)^(ra,*a)) (8) 

and the anomalous 

(n, h; r 2 , i 2 ) = -i {T^{v x ,t x )^t' ( r 2, t 2 )> , (9) 

Green's functions pi Il5j|. 

If the Hamiltonian is time independent, which is the 
case we want to discuss, the Green's functions depend 
on the combination t\ — i 2 not separately on t\ and t 2 
(i.e., G a ^{Y 1 ,t 1 \Y 2 .t 2 ) = G a . a '(r 1 ,r 2 ,t 1 -t 2 ), and sim- 
ilarly for F). Correlation functions (U]) and © can be 
calculated from G and F by the limiting procedures 



/V,<r(r',r) = — i lira G CT / )CT (r',r,£ 



(10) 



Xa',a(r',r) = -tlimi?* (r'.r.-e). (11) 

£->0 

We consider the problem of singlet Cooper paring. In 
that case the nonvanishing elements of the Green's func- 



tions can be chosen to be G-|-|- = G±± and F- 



n 



■it 



respectively [191 ] . For practical purposes let us introduce 
the functions 

u(r 1 ,r a ) = (-— V? + I7(n)-/i U( ri -r 2 ), (12) 
\ 2m ) 

A CTCT '(ri,r 2 ) = v{v x - r 2 )x CTCT '(ri,r 2 ). (13) 
Then, the time evolutions of the two Green's functions 
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can be written as 
8 



ih—G n (ri,r 2l h - t 2 ) = M(ti - t 2 )S(r 1 - r 2 )- 
+ J d 3 rv(r 1 ,r)G n (r,r 2 ,t 1 -t 2 ) 
+ J d 3 rA n {r 1 ,r)F n {r,r 2 ,t 1 -t 2 ), 



(14) 



and 



^^-•P!|.t( r i' r 2)*i - *a) = 
[ d 3 rAl t (r 1 ,r)G n (r,r 2 ,t 1 -t 2 ) 
-J d 3 riy(r 1 ,r)F^(r,T 2 ,t 1 -t 2 ). 



(15) 



The symbol '*' denotes complex conjugation. Let us take 
the Fourier transform with respect to time of the Green's 
functions as 



/oo 
dte wt G n {v u v 2l t) 
-OO 



(16) 



(and similarly for F). Next we transform quantities 
like A(ri,r 2 ) in Eqs. (fT4"]) and (fTSj) to mixed position- 
momentum representation by introducing R = (ri + 
r 2 )/2 and r = ri — r 2 and taking the Fourier transform 
with respect to r: 

A(R, p) = J d 3 k e ipr/h A (R + r/2, R — r/2) . (17) 

We use the term phase space for the (R, p) space in the 
following. If a quantity C(ri,r2) is given by 



C(ri,r 2 ) = J d 3 rA(r u r)B(r,r 2 ) 



(18) 



then Baraff and Borowitz [3 H3] showed that the cor- 
responding relation in position-momentum space can be 
expressed as 



C(R,p) =e[A(R,p),£?(R,p) 



(19) 



where is a bilinear operator acting on two phase space 
functions as [2l| 

G[A(R,p), J B(R,p)] = lim 

R R 



where G = G n (R, p, u), F = F it (R,p,uj), A = 
A|j_(R, p), v = u(R, p), respectively. In deriving Eqs. 
([21~T) . (f22j) we used the properties 

A CTff /(ri,r 2 ) = -A CT / CT (r 2 ,ri) 
A it (R,-p) = A n (R,p) 

which can be proven from the definition of A (Eq. (|13[) ). 
The lu independent functions v and A in the mixed rep- 
resentation are 



i/(R,p) = |- + U(R) - /x 
Zm 



(23) 



and 



A(R,p) = /Ae-.W VrMR + r/2 , R _ r/2) . {2i) 

Green's function are useful for calculating physical 
quantities such as the density n(R) = n-f-(R), or the 
equal-time expectation values /i(R, p) = h^f(R, p) and 
X(R, p) = xn( R 7P) (defined in Eqs. © and flU). From 
Eqs. ([TO]) and {HJ follow that 



n(R) 



d 3 p 
(2tt?i) : 



■/i(R,p), 



fc(R,p) = -« / — G(K,p,w)e ius , 
X (B.,p)* = -i J ^-F(R, p, u)e iue , 



(25) 



(26) 



(27) 



where £ is an infinitesimally small positive regularization 
parameter. 

The widely used interaction potential (JSJ leads to di- 
vergence in the gap equation, which requires some special 
care. Due to the S interaction A(R, p) is momentum in- 
dependent. In that case, the self-consistent equation for 
the local gap A(R) has to be regularized. It means, that 
we should take the regularized part F reg of F by which 
the the self consistent gap-equation 



A(R)* = 



4tt^ 2 c 



provides a finite value for A(R) |22t l23j. 



exp 



ih >— ^ f d d 



d d 



A(R,p)B(R',p'). 

(20) 



Eqs. (fT4j) and (jT5() in the R, p,w representation can be 
written in the compact forms 

h = TiujG — Q[v, G] — 6[A, F], (21) 
Q = huF + Q[v,F]-G[A*,G\, (22) 



III. PERTURBATION SERIES FOR G AND F 

In this section we shall construct a formal solution of 
Eqs. ([21]) and ([22]) supposing that the functions v and 
A are known. The bilinear operator as defined in Eq. 
(I2U1) can be expanded as a formal series of %: 



Q[K 1 ,K 2 ]=Y J ti®j[Ki,K 2 ] 

3=0 



(29) 
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The first two operators 0o and 81 are simply 
9 [K U K 2 ] = K 1 -K 2 , 

Q 1 [K 1 ,K 2 ] = ^{K 1 ,K 2 }, (30) 

where {. . .} is a usual Poisson bracket. For higher order 
terms in the series (|29|) it is useful to treat derivatives 
according to the phase space variables on equal footing 
by the definition 

( d d d d d d \ 

^^^{M^'M^'dR^'dfi'dp^'dp^)' (31) 



where 



We also need an antisymmetric metric g a ° , where 5 14 = 
g 25 = 5 36 = - g 41 = - 5 52 = - g 63 = 1 and all the 
other elements are zero. The metrics reflect the simplec- 
tic structure of the phase space. For example 

6 2 [Ki,K 2 ] = -\g afi g~< s {d^Kx) {d p d 5 K 2 ) (32) 

o 

Expressions for higher order 0j's can be derived simi- 
larly in a straightforward manner. Let us write now the 
normal and the anomalous Green's functions G and F as 
a formal power series in % 

00 

G(R,p,w) = Rj^«, i G i (R,p,w) (33) 
3=0 



F(R,p,o;) = /i^^> i (R,p, 

3=0 



If we write 



(34) 



(35) 



and treat this quantity as an o(h°) term then we get 
from (|2l"j) and (|22|) in different orders of % the following 
equations 



(f2 - v)Gj - AFj = Qj, 
-A*G j + {Sl + v)Fj = Pj, 



with 



and for j ' > 1 



Qo = l, 

Pa = 0, 



(36) 
(37) 

(38) 
(39) 



Qj = E ( 6fc k G ^- fc ] + Gk t A ' ) ( 4 °) 

fe=l 

j 

p l = E t A *' G J-fc] - e * ["1 F J-fc] ) • ( 41 ) 



fc=l 



It is clear from this structure that Qj and Pj for fixed j 
are given in terms of lower order corrections of G and P. 
Solutions to (|3"61 and (|3T|) are 

1 



f2 2 -£ 2 
1 

n 2 -E 2 



[(A + + AP,-] (42) 
[(n - + A*Qj] , (43) 



£ = £(R,p) = vMR,P) + |A(R, P )| 2 



(44) 



Using Equations (|38[) - (1441) one can calculate corrections 
to G and F up to arbitrary large orders. 

Up to now, we have not addressed the question of the 
correct pole structure of G and F. This requires to intro- 
duce infinitesimal imaginary parts in the denominators of 
G and F. This step can be easily performed if we write 
the corrections as partial fractions in fi with 17 indepen- 
dent numerators and choose iS accordingly to 



F; 



k 

E 



A jik (R, p) 
(hu-E + iS) k 

Ci,fc(R,P) 



B i>fc (R,p) 



(huj + E-i5) k 

^i.fc(Rp) 



_(huj~E + iS) k 
The zeroth order coefficients are 

^0,i = ^(l + ^), B 0A - 

and 

Cqa = —Do 1 = 



(hu + E-i6y 



2 V EJ' 



A* 
2E 



(45) 
(46) 

(47) 
(48) 



All the other coefficients are zero. Non-vanishing first 
order corrections involve Poisson brackets in the combi- 
nations of 



~8£3 



(A*{v, A} - A{v, A*} + u{A, A*}) 



(49) 



Ai, 2 - EBi t i+—{A, A*}, B li2 = EB hl -— {A, A*}. 

n ^ (50) 

See also Ref. It is important to note that for real 

A the first order correction G\ to the normal Green's 

function is identically zero. Coefficients of F\ are nonzero 

even for real A as can be seen from 



4£ 



{i/, A*}. 



There are no first order poles of (|46j) for j 
quently, 

d,i =x> M = 0. 



(51) 



1, conse- 



(52) 



Higher than first order coefficients require tedious cal- 
culations. Here we do not give explicitely the sec- 
ond order coefficient functions in the numerators of 
Eqs. (1431) . (|4"rJ)) . Instead, we sketch the structure of 
these corrections. G 2 and F 2 involve k = 1,...,4 and 
the coefficient functions for real A are linear combi- 
nations of ten (usual and) generalized Poisson brack- 
ets {A; v, v}, {A,A} + , {A; A, A}, {A;^,A}, {^;A,A}, 
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{^,A} 2 , {v : v} + , {v;v,v}, {v;A,v}, {v,A}+. The first 
generalized Poisson bracket is denned as 



{A,B}+ = g a Pgi s (d a d 7 A)(dpd s B) 



(53) 



and is symmetric if one makes the changement Ah JS. 
The second generalized Poisson bracket acts on three 
phase space quantities as 

{A; B, C} = sf x(i 9'\d a d 1 A)(dfiB){d 6 0). (54) 



IV. GRADIENT EXPANSION OF PHYSICAL 
QUANTITIES 

In the previous section we have seen that the one parti- 
cle Green's function G can be written as a formal power 
series in h, where the correction terms Gj in (|33l) are 
given by the partial fraction series (|45p . Performing the 
w integal in Eq. ([25)1 it is easy to see that only the first or- 
der poles located on the upper half of the complex omega 
plane give contributions to ft,(R, p). Correspondingly, by 
Eq. ([25]) the density n(R) has the expansion 



n(R) 



OO 

E 

3=0 



d 3 p 
(2nh) 



£jm(R, P) 



3^3 



5>^(R) (55) 

3=0 



Similarly, Eqs. (|34|) . (|46|) and ([27]) and ([24]) lead to 
A(R.p) E^^,(R,p) 

3=0 



3=0 



(56) 



Calculating the first few gj(R)'s and /j(R, p)'s it can be 
seen that the R dependence enters in gj and fj through 
the quantities U(R) and A(R, p) and through the spa- 
tial derivatives of order < j of U(R) and A(R, p). For 
j = there are no spatial derivatives (See Eqs. ([47]) 
and ([48])). For j = 1 the Poisson-brackets in ([49]) bring 
the dependence also on gradients of U(R) and A(R, p) 
into B\ i, and correspondingly into <7i(R) for complex A. 
Bi.i vanishes if A real, and we consider in the following 
only this case. Here we write the j ' = 2 results expressed 
in terms of the generalized Poisson brackets 

„ , 2^ 3 A-3i/A 3 r . 3A 2 i/ r . 
B 2 ,i(R, P) = TTTjyr i A ' ^ ~ ^TE* *4+ 



16.B 7 

2^A 3 -3^A 2 A1 3z^ 2 A 2 
-{A; A, A} 



32S 5 



16£ 7 
A 4 - 4z/ 2 A 2 



8£ 7 



-{A;^,A} 



_ v{1v 2 - 3A 2 )A , . 

w ^ ^7 w A > v i 



8E 7 



16E 7 

5v 2 A\ A1 (2^ 2 -A 2 )A r 

^ A > A ) + Tgcis A )+ 



16£ 7 



16£ 7 



>,A} 2 - 



16£ 5 
2^A 2 



32£ 5 



{ A , A 1^ 



(57) 



(3A 3 - 2^ 2 VA r , 3i/ 2 A rA A1 



16E 7 



A 2 -2z/ 2 A r . ^2A 2 -3^ 2 A rA _ 

■- ^— {i/, z/}+ = — {A; A, ;y} 

32E 5 8E 7 x ' ' s 

v{2A 2 -Zv 2 )A K^ 2 -2A 2 ) r 

-{i/; A, A} — {i/, A}+ 



IQE 7 
5u 2 A 2 



{A; is, v) 



IQE 7 

v i _ 4l/ 2 A 2 



16-E 



— i A ; A , A l 



8E 7 
A 



-{u;A,u} 



16E'- 



>,A} 5 



(58) 



For 5 interaction Eq.([56]) simplyfies, fj(p, R) has no mo- 
mentum dependence: 

OO 

A(R) = X>'£( R ) 



3=0 



AirTi 2 a • °° 



E^/(2W^ (R ' P) ' (59) 



where D - x denotes the regularized part of -Dj,i, which 
can be obtained from Eq. ([56"]) if the pseudo-potential is 
used for the interaction. Equations ([55]) and (|59[) can be 
solved perturbatively whose formal solutions become of 
the form 



i(R) = ^V'nj(R) 



i=o 



A(R) = ^^A,-(R) 

3=0 



(60) 
(61) 



It is important to stress, that on the right hand sides 
of Eqs. ([55]) and ([59]) all the quantities gj and fj depend 
on the total A(R) and U(TL), thus rij ^ gj and Aj ^ fj. 



A. Local density approximation 

The leading order j = approximation in Eqs. ([55]) 
and (|59p are equivalent to the LDA. In that approxima- 
tion one has to solve the equations 



n (R) = <7o(E/o(R),A (R)) 
A (R) = /o([/ (R),Ao(R)). 



(62) 
(63) 



The Hartree-Fock terms in U(R) are density dependent. 
By the notation Uq in that case we mean that they are 
evaluated using no(R). If the Hartree-Fock terms are 
neglected Uq = U ext and Uj = for j > 0. Let us 
introduce the local chemical potential a by 



a(R) = (j,-U(R). 



(64) 



The quantity z/(R, p) defined in Eq. (f2"3"|) is simply 
i/(R, p) = p 2 /2m — a(R). The phase space quantity 
E(R, p) in Eq.(|gi]) with real A(R) is equal to E(R, p) = 



FIG. 1. The dimensionless function ji(t) denned in Eq. (|A4[I . FIG. 2. The dimensionless function ii(t) defined in Eq. (|A3 



s/ ^ 2 (R, p) + A(R) 2 . The go function occuring in ([62 
can be calcultccl from (f55j) and (jT7|) and it is given by 



d 3 p 



r*)» V 1_ ^(RTi)J- (65) 
The momentum integrals can be performed analytically 



^(R),A(R)) = i/^ 



50 (C/(R),A(R)) = 



1 

4^ 



2raa(R) 



3/2 



A(R) 
a(R) 



(66) 

where we have used the dimensionless function ji (x) (See 
Appendix [Al and figH]). The /o function in (|63|) can be 
calculated in a similar way: 



/o(t/(R),A(R)) 



A(R) f d 3 P 



{2nhf 



£(R,p) 



2m 
p2 



The second term in the integrand ensures a finite value 
for the momentum integral, i.e., D i is regularized with 
this term substracted. The momentum integral in Eq. 
(|67p can be written in terms of complete elliptic func- 
tions (see Appendix fA} and can be expressed for negative 
scattering length as 



/o^.ACR^AW^lJ 2 ^-^^ 



h 2 11 V a(R) 



(68) 

where the dimensionless function ii(x), depicted on 
fig. is defined by Eq. (|A3I) . The overall constant chem- 
ical potential /i is fixed by 



N 
~2 



n (R)d 3 R. 



(69) 



Here N is the total particle number (including both hy- 
perfine states). Solutions to Eqs. ([62]), <(63]) and (JSSJ) 
with g and f given by (|6"6"|) and (p5|) are the solutions 
in leading order. Thus, the zeroth order terms in the 



gradient expansion lead to the Local Density Approxi- 
mation. This corresponds to he Thomas-Fermi approach 
generalized to taking into account the pairing field A(r). 
In the following we shall calculate corrections to go and 
/o- 



B. 0(K) order 

We have shown that for real A the quantity Bi i is 
zero. Thus, .gi(R) = in Eq. (|55|) . Due to the property 
(f52|) the regularized part of Di : i is also zero. It means 
there are no corrections to the density and to the gap 
equations in this j = I order. 



C. 0(H 2 ) order 

The evaluations of the j = 2 second order corrections 
<?2(R), /2(R) are rather tedious. In case of momentum 
independent gap and v(R, p) given by Eq. ([23)) nonvan- 
ishing generalized Poisson brackets are 



{.,A} 2 = J2 --(VA),^, (70) 
* — ' mm J 



{u, v } + = -(y 2 U), {^,A} + = 1(V 2 A), (71) 
m m 



m — ' 



d 2 U \ pi pj 



KA, y }=-(VA) (Vt/). 



1 3 
{A;i/,z/} = — V 
m L — ' 

{^;A,A} = 1(VA) 2 . 
m 



d 2 A 

8R~dR, 



HI)', , ) in in 



(72) 



7 



Note that the Laplace operator will be written as V 2 to 
avoid confusion with the gap. We gave a general expres- 
sion for B2.1 in Eq. (|57j) by which g 2 can be obtained 
by evaluating the momentum integrals. For the details 
see Appendix [A] The result is rather lenghty and can be 
presented in the following way. Let us define 

A(a) = A(a(R)) = v 



2TT 2 fl 3 

and the dimensionless combination t by 

A(R) 



t = t(R) 



a(R) 



(73) 



(74) 



where a(R) is the local chemical potential The 
second order correction g 2 to the density equation (|55p is 



+ 



92 (R) = 
(VA) 2 



ma- 



= A(a) 

3 



(V 2 U) 
ma 2 

(V 2 A) 
ma 2 



(VA)(V10 



(Vf7) 



mo. 
2 



-H 5 (t) 



(75) 



Functions H\(t), . . . , H 5 (t) are given in Appendix |B] by 
Eq. ()Blj) . The second order corrections to the gap equa- 
tion can be calculated using Eqs. ([55)1 and (j5"9")l . In sec- 
ond order the momentum integrals exsist, consequently 
D r 2 ® = D 2,1, i.e., there is no need to regularize L>2,i- Pro- 
ceeding as above, the momentum integrals can be treated 
as in Appendix [A] Second order gradient correction f 2 
to the gap equation (I5T)|) can be expressed as 



/ 2 (R) = gA(a) 
(VA) 2 



(v 2 m r/ , (VA)(vto l/M 



TnCt 



ma 

. M 3 (f)+^A/ 4 (t)+^M 5 (() 

mOi mnn,3 



ma° 



(76) 



where the functions Mi(i), . . . , Ms(i) also can be found 
in Appendix|B]in Eq. ([H5]) . From Eqs. (JT5J) and (JTBJ it is 
obvious that the second order corrections to the density 
and to the gap equation involve the spatial derivatives of 
the external potential U ex t and the gap profile A. If the 
generalized Hartree-Fock approximation is considered g 2 
and f% will contain terms with the spatial derivatives of 
the density, too. 

In second order approximation n no + h 2 n 2 , A ps 
Ao + ?i 2 A2, U fa {U ex t + gn>o) + fi 2 gn2- Expanding both 
sides of Eqs. ([55} and (|55|) up to second order in ?i, the 
zcroth order terms cancel. The second order gradient 
corrections to density and the gap are the solution of the 



1 - di9o{U , A ) -<9 2 5o(^o, A ) V n 2 \ _ f g 2 
-dif (U ,Ao) l-d2MU ,Ao)J\A2j \f 2 

' (77) 

inhomogeneous linear equations (Here 9i and $2 denote 
partial derivatives with respect to no and Ao, respec- 
tively) . 

The spatial derivatives of the density in Eqs. (|75j) and 
(f76|) are however missing if the MF-BCS model is consid- 
ered, which neglects the Hartree-Fock terms in v. The 
density and the gap still get gradient corrections in this 
model, which case will be studied next. 



V. SECOND ORDER CORRECTIONS 
CALCULATED PERTURB ATIVELY IN THE 
MF-BCS MODEL 



In MF-BCS model the quantity 

a(R) =H- U ext (R) 



(78) 



is density independent and it is advantageous to use t de- 
fined in (|74"j) instead of A. We keep the density equation 
(1551). but rewrite (f59j) as 



(79) 



3=0 



fi vanishes as in the previously discussed general case. 
Similarly to (fBTj) we are seeking the solution for t(R) as 
a formal series in Ti 



t(R) = 

where t\ is zero due /1 = 0. The first /o is given by 



(80) 



fa = ~\a 



2 /2ma(R)\ 



1/2 



10 1 o 
7T 1 V ^ 



and / 2 can be obtained from Eq. (I76[) by dividing both 
sides by a. 50 is still given by (l66l) 

In the MF-BCS model the leading order LDA equa- 
tions are 



1 /2ma(R) 



no(R) -4^ 2 V ^ 2 

2 /2ma(R) 



3/2 



ji(*o), (82) 



1 = 

7T 



1/2 



ti(*o). (83) 



For fixed chemical potential u the fo(R) profile can be 
calculated from ([53"|). 

To obtain the second order gradient corrections n 2 and 
t 2 we can approximate in the expressions (J75J), (|76p for 
172 and f 2 the quantity f by its zeroth order value to, 
because g 2 and f 2 are already of second order. Taking the 
gradient of the leading order gap equation ([83]) VA(R) 
can be approximated as 

(VA (R)) - -(VEWR))Ti(* ), (84) 
where the dimensionless function T\ (t) is given by 

h(t) 



Ti(t) 



ti 3 (ty 



(85) 



(See Appendix lAl for the definitions of i n (t) and j n (t)) 
Taking the divergence of the two sides of Eq. 
V 2 Aq(R) can be reduced to 



V 2 A (R) = -(V 2 f/ e:ct )7 1 i(to) 



(VU ext ) 2 



T 2 (t ), (86) 
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where we have introduced an other dimensionless t- 
dependent function T 2 (t). Explicitely: 



T 2 (t) 



3i 2 (j 5 (t)i 3 (t)- j 3 (t)i 5 (t)) 



t 2 -il{t) 

(87) 

Using Eqs. (|84|) and (|86|) in the expressions of the sec- 
ond order corrections (|75l) and (|7d| and the results of 
Appendix lAl and IB1 the perturbatively calculated correc- 
tions are 



92 (R) = 



A(a) 



(Vt/ e 



-Pi(to) + (V 2 i7 e:ct )Qi(to) 



and 



/ 2 (R) 



(88) 



gA(a) 



ma* 



(yu e 



-p 2 (to) + (y 2 u ext )Q 2 {t ) 



(89) 

where A(a) is defined by (|73[) . Explicit expressions for 
Pi, Qi, Pi and Q2 are given in Appendix [Bl in Eqs. (IB3p ~ 
(|B6I) . Corrections (|581 . (I59"| involve terms proportional 
to (VU ex t) 2 and V 2 £/ e xt- We remind the reader that g is 
proportional to the scattering length a (see Eq. ([5])). 
Expanding both sides of Eq. (fT9")l up to second order 

using / (t) « /o(i + fi 2 i 2 ) « /o(to) + h 2 t 2 f^(t ) the 
zeroth order terms cancel, and i 2 can be expressed as 



i 2 = - 



1 



M'i(*o) 



2ja| ^ 



2ma(R) 



1/2 



(90) 



Note that t 2 depends on the scattering length only 
through to- The second order gradient correction of the 
gap is A 2 (R) = i 2 (R)a(R) in the MF-BCS model. Us- 
ing the same approximation for the density the second 
order gradient correction of the density is 



n 2 



1 / 2toq;(R) 
i^ 2 I H 2 



3/2 



hj[(to 



92 



(91) 



These are the first nontrivial gradient expansion terms. 
The simplification in the MF-BCS model has arised from 
the fact that in Eq. (f77f the d\ derivatives (i.e., the 
derivatives with respect to the density) are zero. 



VI. UNIVERSAL PREFACTOR OF THE VON 
WEIZSACKER TYPE CORRECTION 

The present paper supply the derivation of some of 
the relations used already in our earlier paper • To 
compare with the results of fiol ] one has to apply the limit 
a — > 00 to Eqs. 



I l91j) (see also Appendix B in applying 
this limit). The leading order to(R) profile is constant, 
which can be seen from Eq. (|83|) . Let us denote by T the 
root of ii(T) = 0, then 



i (R) = T « 1.1622, ji(T) w 1.4688 



(92) 



(T is defined by the requirement ii(T) = in order that 
Eq. (|83[) remains meaningful in the limit a — > 00. See 
also (|B7|-(lB9| . ) It leads to 

2n(R) = ( A1 -C/ e ^) 3/2 2^Ji(T)^) 3/2 

(a) 21 v2 ^« 

2m 3l[ 'ASbi-Uest)* 



'2m Jl( > 192 ~ U ext )3 



(93) 



Note that n = rif = in our notation. Similarly from 
Eqs. GU, (USD, (ESI), USD) one gets 



^(R) — ^(M ~ ^ext) — 



4 + 7r 2 ft 2 v 2 u ext 

36T 2m (ju - U ext ) 



(94) 



7r 2 - 8 h 2 (wu ext ) 2 

144T 2m(/i-f/ e *t) 2 ' 
Eq. (p?3"|) can be rewritten as 

^ Uext = {jl{ T))-y3(2n 2 ) 2 /^(2n) 2 / 3 

7h 2 \ {V{v-U ext )) 2 V 2 {fx~U ext ) 
+ 72m\ 4(n-U ext )z (v-U ext ) 

(95) 

where we have taken into account that the seecond term 
in the right hand side of the equation is a correction. The 
first step of the iteration on the right hand side leads 
after a rearrengement to the Thomas-Fermi- Weizsacker 
type equation: 

7j2 y2 n l/2 

*w - — + £ k f (2n) 2 / 3 + U ext = n (96) 



Here 



K W 



2m n 1 / 2 



27' C 



3ji(T) 



2/3 



Kp 



— (37T 2 ) 2 / 3 . 

2rn ' 



(97) 

Note that £ is the usual universal constant introduced 
for the homogeneous system by the definition /1 = £ep. 
(ep being the Fermi energy h 2 k F /(2m), where kp = 
l&^n) 1 ^). The above value is valid in the MF-BCS 
model @, GjJ. Numerically £ = 0.59 in this model, 
while the Monte Carlo Simulations have provided £ = 
0.37 — 0.44 [2I |24[ . Note, that for a normal system at uni- 
tarity £ = 0.55 (25H27I ]. the corresponding MF-BCS value 
is £ = 1, i.e. free gas value, since £/v = a + /3, where a is 
the ratio of the mass and the effective mass (beeing unity 
in the MF-BCS model) and j3 is zero (see Ref. [27| and 
references therein) . The first term on the left hand side of 
Eq. (j9l)|) is of the form of the von Weizsacker correction 
to the Thomas-Fermi theory (see for the early history of 
the problem Ref. [28[ ) . By now it is well established that 
K w— 1 (originally derived by von Weizsacker) is the cor- 
rect value in case of a rapidly varying density with a small 
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P 

FIG. 3. Solutions of the gap equation without (solid line) 
and with (dashed line) the second order gradient corrections 
as a function of the dimensionless radius p — R/Rtf, and 
measured in units of hu)o = hv (see text). Parameters are: 
a/d = 1/(3tt), R TF /d = 8. 
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P 

FIG. 4. Density profiles without (solid line) and with (dashed 
line) the second order gradient corrections as a function of the 
dimensionless radius p = R/Rtf, and measured in units of 
1/d 3 Parameters are: a/d = 1/(3tt), Rtf /d = 8. 



amplitude, while in case of a smooth external potential 
Kw = 1/9. This value of kw was first derived by Kirzh- 
nits [2!| and by Kompaneets and Pavlovskii [30j . (S ee for 
reviews of the density gradient expansions [ll], |3l|, |32|). 

It is worth mentioning that kw = 1/9 was found [331 ] 
the optimal value when the energy of a free gas in a har- 
monic oscillator potential was compared with the quan- 
tum mechanical result via second order perturbation the- 
ory. This suggests that such an external potential occur- 
ing in trapped gases is well suited for a gradient expan- 
sion of the density. 

There has been a renewed interest in recent years 
concerning the von Weizsackcr correction in case of the 
trapped unitary Fer mi gas |34l443l |. The value of kw = 1 
has been chosen in 1341 l3al ■ while kw = 1/4 has been 
obtained in (3(| H3, EUa, Eo| by assuming the validity of 
a kind of Ginzburg-Landau theory at zero temperature. 
Furthermore, an expansion in powers of d — 4 — e spatial 
dimensions has led to Kw = 0.176 42[ by extrapolating 
the result to three dimensions. A comparison between 
the choices k = 1/9 and ft = 1/4 has been carried out in 
43] by studying fermion systems at unitarity with parti- 
cle numbers up to 50. It has been found that the choice 
kw = 1/4 provides better results for the energy except 
at few particle numbers. This finding backs our result 
for kw = 7/27, which is quite close to this value. 



VII. ISOTROPIC HARMONIC TRAPPING 



As an application of Eqs. (|55| - (|9"Tj) let us apply our 
results to the special case of isotropic harmonic trapping 
potential 



U(R) 



-mujlR 2 . 



(98) 



In local density approximation the Thomas Fermi radius 
Rtf is introduced by the relation 



1 



muinR 2 



(99) 



which ensures uo(Rtf) = 0. It is advantageuos to use 
the dimensionless combination 



g — Rj Rt f 



(100) 



for the radial distance. A natural characteristic length of 
the ha rmonic os cillator problem is the oscillator length 
d = y/h/ (mwo). The LDA gap equation (|83"1) for har- 
monic confinement 



= 2 | d | Rtf 

7T d 2 



(l-e 2 ) 1/2 n(io) 



(101) 



provides us a profile to(g), which depends on the single 
dimensionless parameter \a\RTF/d 2 ■ Up to second order 
gradient corrections the gap can be expressed as 



A(e) 



where 5t can be read off from (|90p as 

st= 1 fWg 2 P 2 (t ) , 24Q 2 (t ) 



*o<i(*o) V (1-f 2 ) 3 (W 2 )' 



(102) 



(103) 



Eq. (|9l"T) together with the leading LDA for the density 
can be expressed as 



d 3 n{g) 



(1 - g 2 ) 3 / 2 R 7 



At: 2 



(P 



J'i(to) 



R 



TF 



(104) 
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where 



Sn(g,t ) 



I6g 2 



+ 



(1 - P 2)3 

24 



toi'itto) 

2,2 »Ql(*0) + 74^rQ2(/ul)J0.-,* 



*o«i(*o) 



It is clearly seen that the small parameter of the prob- 
lem is (I/Rtf- The magnitude of the correction as com- 
pared to the leading term is proportional to (cI/Rtf) 4 
both for the density and the gap. 

At the Feshbach point (a — > oo) our results can be 
further simplify 



St 



Sn 



IT 2 - 8 q 2 
36T (l-£? 2 ) 3 



7T 2 



1 



6T (1 - g 2 ) 2 



-7ji(T) 



1 



24(l-g 2 ) 3 4(l-g 2 



, (106) 



(107) 



St and Sn are universal at unitarity (at the Feshbach 
resonance) for a spherical parabolic trap: they do not 
contain any parameter of the two particle interaction. 



VIII. SUMMARY AND CONCLUSIONS 

We have calculated the gradient corrections on the 
BCS side of the Feshbach resonance to the generalized 
Thomas-Fermi model, which represents the LDA in the 
presence of pairing. Though the correction terms have 
a prefactor, which is small for typical trap potentials al- 
ready at moderately large particle numbers, the correc- 
tions get large due to the singularities at the LDA border 
of the cloud. At unitarity a von Weizsacker type correc- 
tion appears whose universal prefactor has been derived 
as kw — 7/27. This value is quite close to 1/4 pro- 
posed in refs. [36|, l39M4l1 | and is also not far from the 
e-expansion result [42[ as extrapolated to three dimen- 
sions. 

It is remarkable that by inverting the functional n[[/ e;r t] 
to order h 2 , as it has been done in Sec. I VII the singu- 
larities at n — U ex t disappear and the density can be 
continued to infinity. This situation is similar to what 
happens in case of the free gas, and perhaps, the most 
physical justification is, which starts the calculation at fi- 
nite temperature and the zero temperature limit is taken 
at the end [Il|,[44],[45j]. One has to keep in mind, however, 
that it does not mean that even the asymptotic decay of 
the density follow the true one in general. 

Away from unitarity, however, the situation is much 
more complicated and needs further study. Instead then 
one can use the treatment applied in Sec. IVIII (i.e., to 
regard the gradient terms as corrections and keep away 
from the surface region, which becomes, however, larger 
and larger when tending to the BCS limit). 



In Figure [3] we have depicted the gap profile A(g)/hu>o 
both in LDA and with gradient corrections. At a certain 
radius g2 the gap with gradient corrections becomes zero. 

In Figure 2] we show the dimensionless density profile 
at the same parameters as for Figure [3J The deviation 
from the LDA profile is much less pronounced at those 
particular parameters in the region where the gradient 
expansion is applicable. Note that in the figures both 
curves are calculated at the same \i values, so they belong 
to slightly different particle numbers. 

The distance g\ from the origin, where A(gi)/hujQ = 1, 
decreases when the magnitude of the scattering length 
becomes shorter, which means that the most suitable sit- 
uation exists at the Feshbach resonance. In the weak cou- 
pling (BCS) limit A is smaller than huo already at the 
point r = 0. For g > g\ the A(g)/huJo function steeply 
goes to zero (see Fig. [3]) beyond which point even its 
formal continuation becomes meaningless reflecting the 
fact that such an expansion is not adequate when the 
gap function A(r) gets smaller than the level spacing of 
the trap. One has to emphasize that this behavior has 
been shown when the first nonzero correction is treated 
perturbatively. More generally, the solution levels off for 
increasing g and one can define the radius gi in such a 
way A(g)/hu}Q < S for g > g 2 with S as a suitable chosen 
small parameter. Actually, in the region g > g 2 one has 
to apply another method instead the one developed in 
this paper to get more accurate result, but the difference 
might be small. This problem goes beyond the scope of 
the present paper and planned as a forthcoming work. 

Note added: After submission of the paper we have 
learned that the density matrix in case of the inhomoge- 
neous superfluid Fermi systems was derived in Ref. [461 ] 
to o(H 2 ) using the Wigner-Kirkwood ^-expansion method 
by regarding the pairpotential as an external one, which 
is an intermediate step in our work (see also 0). We are 
greatful to prof. Schuck for informing us of the papers 
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Appendix A: Momentum integrals 

In the zeroth order two types of momentum integrals 
occur: 

f d 3 p (\ 2m\ A(a) , .. /A „, 



(2nh) 3 \E p 



L > = 1 J^W ^ ~ = ^^It^A/a, (A2) 



11 



where A(a) has been introduced in (|73|) with a the local 
chemical potential The dimensionless integrals i\{t) 
and j\ (t) are defined as 



h(t) = dx\ —======== 



- 1 



h(t) = / £ 2 (fa; 1 =. 

W Jo I ^ix 2 -l) 2 + t 2 



(A3) 
(A4) 



and ji (t) can be expressed in terms of complete 
elliptic integrals K(k) and E(k) (see Reference |47| ) 



K{k) 
E(k) 



tt/2 



dip 



1 



tt/2 



yl — k 2 sin 2 </p 
c^a/ 1 — k 2 sin 2 <p, 



as 



»i(f) = v 7 ! + t 2 - 2B(fc)] 



t 2 K(fc) 



where the modulus A: is connected to t by 



2E(k) 



VTTT 2 



(A5) 
(A6) 

(A7) 
(A8) 

(A9) 



In the special case t = 0: fc = 1, ii(0) = oo, ji(0) = 
2/3. In higher orders one needs the generalizations of the 
integrals (|Al|) and (|A2|) . For n = 3, 5, . . . let us consider 
the momentum integrals 



L, 



d 3 p 1 A(cv) 



(2tt?i) 3 E n a r ' 



in(t)\t=A/ a , n > 1, 

(A10) 



d 3 p j/ A(a) ... 

^Jn(*)|t=A/a, n > 1. 



(2tt?i) 3 £™ a ,! 

(All) 

Here the new dimensionless integrals i n (t), j n {t) are de- 
fined for odd n as 



oo 



1 



in(t) = I x 2 dx — , 

' Jo W(^-l) 2 +* 2 



x~ 



1 



n > 1, 
(A12) 

n > 1, (A13) 



j n (t) = / a^efo— , 

■>0 (y^ - l) 2 + t 2 J 

For similar integrals written in a different way see Ref. 
[48|). They can obtained analitically from ii(t) and Ji(t) 
using the rules 



*2„+lW = (-l) r 



1 



l-3---(2n-l) \tdt 



1 Pi \ n 

--) ^t) (A14) 




FIG. 5. The Pi(t) function (see Eq. (|B3 



i2»+i(t) = (-1)' 



i 



1 9 



1 • 3- •• (2n- 1) V* dt 



hit), 



(A15) 

which can be easily seen from definitions (|A12[) and (|A13|) 
respectivelly. Useful properties performing the gradient 
expansions are 

j' 1 (t)=H 3 (t), (A16) 
i' n it) = -nt*„+2(<), (A17) 
j' n (t) = -ntj n+2 it), n>l. (A18) 

In calculating explicitly i n (i) and j n (i) for odd n using 
the well known formuli for the derivatives of complete 
elliptic functions (4tJ it turns out that they are linear 
combinations of ii(t) and Ji(i): 

*n(*) = + B n (t) n {t) (A19) 

jn(*) - tf n (i)*i(i) + (A20) 

where the coefficients A n (i), B n it), C n {t) and D n it) are 
rational functions of i. 



Appendix B: Second order coefficients 

Here we enumerate some dimensionless functions used 
in the main text. Functions occuring in Eq. (|75|) are 



48Hj.it) = 8t 2 i 5 it) - \QtH 7 it) -t 2 j 5 (t) - 10i 4 j 7 (t), 

8H 2 it) = 5t 3 j 7 it)-2tj 5 it), 
48H 3 it) = 2i 3 (t) - 17t 2 z 5 (i) + 15t 4 i 7 (i) +2j B (t), 
48# 4 (*) =2h 3 (f)+5f 3 i 5 (f) - 10t 5 « 7 (i) 

-4<j 5 (t) + 10t 3 j 7 (0, 
16# 5 (t) = 4£ 2 i 5 (t) - 5i 4 i 7 (i), 

(Bl) 



12 



-0.018 



-0.019 



"J -0.0195 

a 



-0.0205 





FIG. 6. The Qi{t) function (see Eq. (|B4 



FIG. 7. The P 2 (t) function (see Eq. {B5 



0.25 



and those used in (|7o| are 



48Mi(t) = 10 1 3 j 7 (t) - 4tj 5 (t) + 2ti 3 (t) 
+5t 3 i 5 (t) - 10t 5 i 7 (t), 
8M 2 (t) = i 3 (t) -5t 2 i 5 (t) + 5t 4 i 7 (t), 

48M 3 (t) = 11 tj 5 (t) - 15t 3 j 7 {t)+2ti 5 (t), 

48M 4 (<) = -3j 3 (t) - t 2 j 5 (t) + I0t 4 j 7 (t) 
-I0t 2 i 5 (t) + 10t\{t), 

16M 5 (t) = 5t 3 j 7 {t)-2tj 5 (t). 



(B2) 




0.05 



FIG. 8. The Q 2 {t) function (see Eq. (|B6 



Straightforward, but lengthy calculation leads to the 
analytic forms of the coefficient functions Pi(t), Qi(t) oc- 
curing first in Eq. 



[8 + 3t 2 ] h{t) 



5ji(t) 



384(1 + t 2 ) 128(1 + t 2 ) 
t 4 [l + t 2 ] ij(t) t 2 i\{t)[t 2 + Z\ 



192 (3 h (t) - P h{t)f 192 (3ii (t) -t 2 *i(t)) 
i 2 {t) [4 + 3t 2 ] 
384 {Z n {t) - t 2 n{t)y 



(B3) 



P*(t),Qi(t) (used inEq. 

nit) 5t[*i(«) + 3ji(t)] 



48 i 384(1 + t 2 ) 

*ii(t)[5*i(t) + 7ji(t)] 



384(< 2 zi(i) -3ji(i)) 
| tjift) [-8 ? 2 ft) +21j 2 (t)] 

384 (t 2 ii(t)-3ji(t)) 2 
+ 3h(t)} t ll (t)j 1 (t) [ lOhjt) + 21j 1 (t)] 
192(t 2 ii(t)-3ji(t)) 3 



(B5) 



W ' ~ 72t + 96(1 + t 2 ) 

Ql (t)= tH{t) h{t) + ^ *ii(t)[10ii(t) + 21ii(t)] , , 

U 32(l + i 2 ) + 96(3, lW -t 2 il W) ) - 2 y^ i( V_3 Ji( V/ J - (B6) 

The functions Pi(t), Qi(t), ^(i) and are shown in 

Similar calculation gives the expressions for Figs. [31 El [3 and [5] respectively. At the Feshbach reso- 



13 



nance ii(T) = should be taken. In this case 



Pi(T) = 



h (T) , 5T n (T) 



48T 128(1 + T 2 )' 



(B8) 



P X (T) 



128(1 + T 2 )' 



Qi(T) 



32(1 + T 2 )' 



(B7) 



Qa(T) 



ji(T) , Tj\(T) 



24T 32(1 + T 2 )' 



(B9) 
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